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ROUTINE - NOPMAN 


PURPOSE - The main program which begins the calculation. 
It also does the summation over radius. 

AUTHOR - Roy K. Amiet 

INPUT 


COMMON BLOCKS 


COMMON 

/ ROTCA/ 


Name 

Type 

Description 

CL 

RD 

Chord/turbulence length scale 

BN 

RD 

Number of blades 

COC 

RD 

Sound speed / chord 

RC 

RD 

Blade radius /chord 

R1C 

RD 

Far-field distance / chord 

RPS 

RD 

Rotational speed revolutions / sec 

DXDZ 

RD 

Deformation tensor 

ALPHA 

RD 

Tip path plane angle of attack 

AO 

. RD 

Speed of sound 

TMIN 

RD 

Minumum observer polar angle in degrees 
(theta=0 along rotor axis in thrust direction) 

TMAX 

RD 

Maximum observer polor angle in degrees 

D 

RD 

Diameter of rotor 

DELT 

RD 

Desired increment in polar angle between 
successive calculations 

PMIN 

RD 

Minimum observer azimuthal angle in degrees 
(phi=0 in direction of mean wind) 

PMAX 

RD 

Maximum observer azimuthal angle in degrees 

DELP 

RD 

Desired increment in azimuthal angle between 
successive calculations 

ZM 

RD 

Axial Mach number 

FM 

RD 

Flow Mach numbr in rotor plane 

UU 

RD 

Rms turbulence velocity / axial velocity 

CHS 

RD 

Geostrophic wind speed (m/s) 

WUINF 

RD 

Vertical component of the rms turbulence 
normalized by the stream velocity 

COMMON 

/NOPCB/ 


Name 

Type 

Description 

J 

I 

Number of azimuthal integration points 

HF 

I 

First harmonic 

HL 

I 

Last harmonic 

NUM 

I 

Number of points 

IH 

I 

Homogenous or nonhomogenous indicator 

0 homogenous case 

1 nonhomogenous 

NR 

I 

Number of radial points 



LOCAL VARIABLES 


* 


* 

* 


* 

* 

Kane 

Type 

Description 

* 

DH 

RD 

Fractional harmonic spacing between frequency 

* 



calculations 

★ 

DR 

RD 

Increment between radial points 

* 

H 

RD 

Frequency of sound in multiples of blade passage 

t 



harmonics 

it 

HOLD 

RD 

Last harmonic calculated 

it 

IARRAY 

I 

Array containing integer words from a data 

it 



member record 

it 

12 

I 

Do loop index 

it 

II 

I 

Do loop index 

* 

IP 

I 

Number of P values 

t 

ISTAT 

I 

Status of calls to data member routines 

it 

JJ 

I 

Do loop index 

it 

K 

I 

Do loop index 

it 

LPOUT 

I 

Logical unit to write 

t 

N0PIH2 

RD 

Array containing data member name of HOP 

* 



input file 

t 

NP 

I 

Do loop index 

* 

HR 

T 

l 

Humber of radial points 

* 

NT 

T 

4 

Do loop index 

it 

HURDS 

I 

Humber of words read in from data member record 

it 

P2 

RD 

Pressure averaged around azimuth 

* 

P22 

RD 

Previous P2 value 

* 

PH 

RD 

Angle of flight Mach number 

t 

RARRAY 

RD 

Array containing the real words from a data memb 

• 



record 

* 

RCP 

RD 

Radial point 

t 

RCP1 

RD 

Previous radial point 

t 

SUM 

RD 

Sum of pressure squared contributions to noise 

* 

t 

TH 

RD 

Polar angle of observer 

t 

£ 

OUTPUT 



* 

* 

Name 

Type 

Description 

* 

H 

RD 

Frequency of sound in multiples of blade 

t 



passage harmonic 

* 

FI 

RD 

Frequency of sound in Hz 

* 

PSD 

RD 

Spectrum level in dB (per unit Hz) 

t 

NP1 

RD 

Summation counter 

t 

NM1 

RD 

Summation counter 

* 

£ 

N 

RD 

Summation counter 

* 

FUNCTIONS 



* 

1. To call 

subprogram to read in input values 

* 

2. To call 

subprogram giving spectrum Integrated over azimuth 

* 

* 

3. To perf 

orm integration of spectrum over radius 

* 

SUBPROGRAMS 

CALLED 


* 

MMCLOS MMGETR MMOPRD 

NOPINP PX XFETCH XSTORE 
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* 

* 

* 

* 

* 

t 

* 

* 

* 

* 

* 

* 

* 

* 

it 

± 

it 

it 

it 

* 

* 

t 

it 

t 

• 

* 

it 

* 

* 

* 

* 


CALLING PROGRAMS 
ANOP EXEC 

ERRORS 

NON-FATAL 

1. Error opening data member for input 
FATAL 
none 

ENTRY 

Print the input quantities for checking 

Write the values of the deformation tensor 

If beginning harmonic is negative, return to matrix input 

If harmonic spacing is negative, return to previous input 

If number of steps is negative, return to first input 

Calculate frequency step size 

Print headings for the eventual output 

Initialize the frequency variable 

Input of zero for L2 defaults to a single frequency step 
DO 60 K=1,L2 

Increment frequency 
Initialize radius 
Initialize 3um 
DO 22 12=1,19 

Increment blade radial position 
Calculate tip Mach number at this radius 

Call subroutine giving azimuthally averaged spectrum at this 
radius 

Halve first value in trapezoidal integral over radius 
Add to integrated spectrum 
END DO 

Calculate PSD 
Calculate harmonic 
Write results 
END DO 
EXIT 
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ROUTINE - NQPINP 


* 

* PURPOSE - To read in input variables and calculate the mean 

* and turbulence properties for an atmospheric boundary 

* layer 

* 

* AUTHOR - J. C. Simonich 

* 

* INPUT 

* 


* USER PARAMETERS 

* 


ft 

* 

Name 


Type 

Description 

* 

AO 


RD 

Speed of sound (m/s) 

ft 

D 


RD 

Diameter of rotor (m) 

* 

C 


RD 

Blade chord (m) 

ft 

BN 


RD 

Number of blades 

ft 

FF 


RD 

Distance from rotor to observer (m) 

ft 

RPS 


RD 

Rotational speed in rev/sec 

ft 

HF 


I 

First harmonic 

ft 

HL 


I 

Last harmonic 

* 

NUM 


I 

Number of frequency points 

ft 

TMIN 


RD 

Minimum observer polar angle in degrees 

# 




(theta=0 along rotor axis in thrust direction) 

* 

TMAX 


RD 

Maximum observer polar angle in degrees 

* 

DELT 


RD 

Desired increment in polar angle between success 

ft 




calculations 

* 

PMIN 


RD 

Minimum observer azimuthal angle in degrees 

t 




(phi=0 in direction of mean wind) 

* 

PMAX 


RD 

Maximum observer azimuthal angle in degrees 

* 

DELP 


RD 

Desired increment in azimuthal angle betveen 

* 

ft 




successive calculations 

* 

ft 

DATA MEMBER 

ROTNOP(ABLOTl) - This i3 the output from the FM ABL 

* 

ft 

GWS , 

LWX, 

WUINF 


* 

ft 

Name 


Type 

Description 

* 

GWS 


RD 

Geostrophic wind speed (m/s) 

* 

LWX 


RD 

Correlation length scale 

ft 

WUINF 


RD 

Vertical component of the rms turbulence 

* 

t 




normalized by the velocity 

* 

ft 

DATA MEMBER 

R0TN0P (R0T0T1 ) - This is the output from the FM ROT 

* 

TITLE 




* 

NR, J 

. IH 



* 

R2, ALPHA, 

R 


ft 

STRENG 



• 

If IH 

=0 (homogenous case) read in following: 

ft 

T, U, 

V, W 

, X, Y, 

2 

• 

XSL1, 

YSL1 

. ZSL1, 

XSL2, YSL2, ZSL2 

ft 

USAV, 

VSAV , WSAV 
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* 

XA11, 

XA21, XA31 


* 

XA12, 

Xa22 , XA32 


* 

XA13, 

XA23, XA33 


A 

* 

ft 

Name 

Type 

Description 

k 

TITLE 

A 

Title of case run 

k 

NR 

i 

Number of radial points 

k 

J 

i 

Number of azimuthal points 

k 

IH 

i 

Indicator of homogenous or nonhomogenous case 

* 

R2 

RD 

Maximum radius 

* 

ALPHA 

RD 

Rotor tip path plane angle of attack 

* 

R 

RD 

Radius of rotor 

k 

STRENG 

RD 

Vortex circulation strength 

k 

For homogenous case (IH=0) the following input is read, otherwise 

k 

the rest of the 

data member is read in the subroutine NOPPX 

k 

T 

RD 

Streamline drift time 

k 

U 

RD 

U velocity component in radii/sec 

k 

V 

RD 

V velocity component in radii/sec 

k 

W 

RD 

W velocity component in radii/sec 

k 

X 

RD 

Downstream X coordinate of streamline in 

k 



rotor plane 

k 

Y 

RD 

Downstream Y coordinate of streamline in 

k 



rotor plane 

k 

Z 

RD 

Downstream Z coordinate of streamline in 

k 



rotor plane 

k 

XSL1 

RD 

Downstream X coordinate of streamline in 

k 



standard coordinate system 

k 

YSL1 

RD 

Downstream Y coordinate of streamline in 

k 



standard coordinate system 

k 

ZSL1 

RD 

Downstream Z coordinate of streamline in 

k 



standard coordinate system 

k 

XSL2 

RD 

Upstream X coordinate of streamline in 

k 



standard coordinate system 

k 

YSL2 

RD 

Upstream Y coordinate of streamline in 

k 



standard coordinate system 

k 

ZSL2 

RD 

Upstream Z coordinate of streamline in 

• 



standard coordinate system 

* 

USAV 

RD 

Uptsream U velocity component in radii/sec 

k 

VSA V 

RD 

Upstream V velocity component in radii/sec 

k 

ZSAV 

RD 

Upstream Z velocity component in radii/sec 

k 

ft 

XA12, . 

RD 

Deformation tensor 

k 

ft 

LOCAL VARIABLES 


t 

t 

Name 

Type 

Description 

* 

IARRAY 

I 

Array to store a data member record 

k 

IHDR 

I 

Array containing the length of the largest 

ft 



record and the number of records written 

k 

II 

I 

Do loop index 

k 

ISTAT 

I 

Status of calls to data member routines 

k 

ITYPE 

T 

A 

Indicator of user parameter type 

* 

JJ 

I 

Do loop index 

* 

LPOUT 

I 

Logical unit to write 
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MEXIST 

T 

X 

Indicator of existance of data member 

NEL 

I 

Indicator of number of elements in user 
parameter 

NOPIN 1 

RD 

Array containing data member name of ABL output 
file 

N0PIN2 

RD 

Array containing data member name of ROT output 
file 

NURDS 

I 

Number of words read in from data member record 

PI 

RD 

3.14159. . . 

RARRAY 

RD 

Storage area for real variables being read in fr 
data member record 


OUTPUT 

COMMON BLOCKS 

NOPCA - sea NOPMAN 
HOPCB - see NOPMAN 
FUNCTIONS 

1. To input variables for NOP. 

2. To input variables from FM HOT. 

3. For homogenous case read in deformation tensor and calculate axial 
mach number, flow mach number and rms turbulence velocity. 

SUBPROGRAMS CALLED 

MMCLOS, MMGETR , MMOPRD , MMVUM , NEUTRAL, STABLE, UNSTAB, XFETCH, 
XSTORE, XGETP , XASKP 

CALLING SUBPROGRAM 
NOPMAN 

ERRORS 

NON-FATAL 

1. Error finding user parameter in table 

2. Error in opening data member 

3. Error reading in record from data member 
FATAL 

none 

ENTRY 

Input values from data member ROTNOP (N0PIN1 ) 

If L > 0 call STABLE 
If L < .0 call UNSTAB 
If L = 0 call NEUTRL 

Calculate vertical component of the rms turbulence normalized by the 
stream velocity 

Input data member ROTNOP (R0T0T1) 

If IH=0 read in deformation tensor 
EXIT 



* ROUTINE - NOPPX 

t 

* PURPOSE - To take the cross product of column vectors N2 and N3 

* of E and place the result in column Ml. 

* 


* 

£ 

AUTHOR - Roy X. 

Amiet 


* 

INPUT 



* 

ARGUMENTS 



* 

'ft 

Name 

Type 

Description 

* 

NQPIN2 

RD 

Array containing data member name 

* 

FI 

RD 

Observer frequency 

t 

RCP 

RD 

Radial point 

* 

TH 

RD 

Polar angle of observer 

* 

PH 

RD 

Angle of flight Mach number 

£ 

COMMON BLOCK 



* 

NOPCA 



it 

* 

£ 

Name 

Type 

Description 

* 

CL 

RD 

Chord/turbulence integral length 3cale 

* 

BN 

RD 

Blade number 

* 

COC 

RD 

Sound speed/chord 

* 

R1C 

RD 

Far-field distance/chord 

* 

RPS 

RD 

Revolutions /second 

* 

DXDZ 

RD 

Deformation tensor 

* 

ZM 

RD 

Axial Mach number 

it 

FM 

RD 

Flight Mach number in rotor plane 

it 

X 

UU 

RD 

RMS turbulence velocity/axial velocity 

M 

* 

COMMON BLOCK 



* 

NOPCB 



* 

* 

* 

Name 

Type 

Description 

* 

J 

I 

Number of azimuthal integration points 

• 

IH 

I 

Specifies homogenous or 

it 

* 



inhomogenous calculation 

* 

OUTPUT 



X 

* 

X 

ARGUMENT 



* 

* 

Name 

Type 

Description 

* 

P2 

RD 

Pressure averaged around azimuth 

* 

NP1 

I 

Counter for upward summation in NTITRB 

* 



subroutine 

* 

NM1 

I 

Counter for downward summation in NTITRB 
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subroutine 

N I Counter for rescaling of CVT in MTITRB 

subroutine 

LOCAL VARIABLES 


Name 

Type 

Description 

TM 

RD 

Local tip Mach number 

AJ 

RD 

Float J 

DEL 

RD 

Azimuthal step size 

C 

RD 

Cosine theta 

S 

RD 

Sine theta 

C4 

RD 

Co3ine psi 

S4 

RD 

Sine psi 

QM 

RD 

Expression in analysis 

BFZ2 

RD 

Expression in analysis 

SQ1 

RD 

Expression in analysis 

RER 

RD 

Retarded radius/actual radius 

SUM 

RD 

Sum of azimuthal spectral contributions 

G 

RD 

Azimuthal angle gamma 

I 

I 

DO loop counter 

C2 

RD 

Cosine gamma 

S2 

RD 

Sine gamma 

C5 

RD 

Cosine (gamma + psi) 

S5 

RD 

Sine (gamma + psi) 

ALP 

RD 

Angle alpha 

Cl 

RD 

Cosine alpha 

SI 

RD 

Sine alpha 

RM 

RD 

Mach number component along chord 

C3 

RD 

Cosine of angle phi in analysis 

X 

RD 

Observer coordinate 

Y 

RD 

Observer coordinate 

Z 

RD 

Observer coordinate 

SG 

RD 

Modified radius 

FP 

RD 

Frequency measured on blade 

XX 

RD 

Value of x wavenumber 

YK 

RD 

Value of y wavenumber 

T 

RD 

Blade passage time 

T1 

RD 

Time between eddy chops 

XX 

RD 

Specific value of x in analysis 

YY 

RD 

Specific value of y in analysis 

zz 

RD 

Specific value of z in analysis 

T2 

RD 

Time T1 plus propagation time difference 

CVT 

RD 

Step size for summation over wavenumber 

ZKO 

RD 

Initial value of wavenumber in summation 

AD 

RD 

Contribution to sum of particular azimuthal 
station 

SUBPROGARAMS 

CALLED 



NOPNI NOPLFT 


CALLING PROGRAMS 
NOPMAN 


ERRORS 



* NON-FATAL 

* 1. Error reading in data member record 

* fatal 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

t 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

it 

* 

* 

* 

* 


none 

ENTRY 

Calculate TM, the local tip Mach number 
Float J 

Calculate DEL, the azimuthal step size 

C, S and C 4, S4, are the cosine and sine of theta and psi 
respectively 

QM is a factor appearing in the analysis 

BFZ2 is a Prandtl-Glauert factor 

SQ1 is a factor appearing in the analysis 

RER is the observer retarded distance/actual distance 

Initialize NP1 , . . . .G 

DO 50 1 = 1, J Integration aver azimuthal angle gamma 

Increment gamma 

Find cosine and sine of gamma and gamma * psi respectively 
Calculate ALP, the angle of the rotor blade vrt the rotor plane 
Calculate Cl and SI, the cosine and sine of alpha 
Calculate RM, the Mach number of rotor segment vrt to fluid 
Calculate C3, the cosine of an angle in the analysis 
Calculate observer coordinates in rotor fixed coordinates 
Calculate 3igma, a modified observer distance 
. Calculate FP, the frequency on the blade 
Calculate XX and YK , the x and y wavenumbers 
Calculate T, the time between blade passes 
Calculate Tl, the time between eddy intersections 
Calculate XX, YY, ZS which are X, Y, Z values 
Calculate CVT which is 2*Pi/eddy passing distance Z 
Calculate ZKO, the initial radian frequency in the summation 
Call turbulence summation subroutine 
Call airfoil response subroutine 
Increment counters 

Contribution to spectrum from azimuthal integration 
Summation over azimuthal spectrum 
End DO 

Multiply by remaining factors 
EXIT 
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* * * 


ROUTINE - HOPCRS 


* 

* 

t 

* 

t 

* 

* 

* 

it 

* 

* 

t 

* 

* 

* 

* 

* 

* 

* 

* 

t 

* 

* 

* 

* 

* 

* 

* 

it 

* 

* 

* 

* 

* 

* 

* 

t 

* 

* 

* 

t 

* 

* 

* 

* 


* 

* 

* 



PURPOSE - To take the cross product of column vectors N2 and N3 of E 
and place the result in column II. 

AUTHOR - Roy K. Amiet 

INPUT 


ARGUMENT 


Name 

Type 

Description 

E 

RD 

3x3 rotation matrix. 

Nl 

I 

Column into which the cross product is to be 
placed. 

OUTPUT 



ARGUMENT 



Name 

Type 

Description 

V 

RD 

Unit vectors along coordinates 

LOCAL VARIABLES 


Name 

Type 

Description 


N2,N3 I Column vectors for which the cross product is found. 

EN1M RD Normalizes the output to a magnitude of 1. 

SUBPROGARAMS CALLED 
None 

CALLING PROGRAMS 
NOPKVC 

ERRORS 

None 

ENTRY 

Calculate N2 and N3 so that Nl, N2, N3 are in cyclic permutation. 
Calculate cross product of column vectors H2 and N3; place result 
column Nl. 

Find the magnitude of vector HI. 

Normalize vector Nl. 

EXIT 
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SOUTINE - HOPKVC 


PURPOSE - To calculate effect of rapid distortion on turbulence 
spectrum. 

AUTHOR - Roy K. Amiet 
INPUT 
ARGUMENT 


Name 

Type 

Description 

DR 

RD 

Wavevector 

DXDZ 

RD 

Deformation tensor 

VN 

RD 

Unit vector normal to blade 

OUTPUT 

ARGUMENT 

Name 

Type 

Description 

ED 

RD 

Three coordinate vectors downstream of 
contraction. 

EU 

RD 

Three coordinate vectors upstream of 
contraction. 

DKM 

RD 

Magnitude of downstream wavevector. 

DKUR 

RD 

Ratio of downstream to upstream wavevectors. 

UQDQ 

RD 

Ratio of upstream to downstream turbulent 



velocities. 

LOCAL VARIABLES 


Name 

Type 

Description 

F(I) 

RD 

For temporary storage of EU(I,3). 

AF1 

RD 

Magnitude of wavevector, to be used for 
normalization. 

FUNCTIONS 

1. To calculate the 

effect of rapid distortion on the wavevector 

and Fourier 

component amplitude 

SUBPROGARAMS 

CALLED 



NOPCRS 

CALLING PROGRAMS 
NOPNI 

ERRORS 

NONE 


ENTRY 

Find magnitude of uavevector. 


BO I = 1,3 

Set third column of ED equal to uavevector direction. 

Set second column of ED equal to blade normal. 

END DO 

Find cross product of first two vectors 

Find cross product of previous result with third vector of ED. 

DO 1=1,3 
DO J=1 ,3 

Initialize EU 
DO K=l,3 

Take product of ED with deformation matrix to find EU. 

END DO 
END DO 
END DO 

Store third column of EU in F. 

Find the magnitude of first column of EU. 

Normalize the first column of EU 

Find the cross-product of columns 1 and 2 of EU. 

Find cross-product of previious result vith third vector of EU. 
Calculate ration of downstream to upstream vavevector magnitude. 
Calculate ratio of upstream to downstream turbulent velocities. 
EXIT 



ROUTINE - HOPFNL 


PURPOSE - To calculate the Fresnel integrals C and S. Ref. Abromovitz 
and Stegun, p 302. 

AUTHOR - Roy X. Amiet 

INPUT 

ARGUMENT 


Name 

Type 

Description 

X 

RD 

Fresnel integral argument 

OUTPUT 



ARGUMENT 



Name 

Type 

Description 

E 

C 

Fresnel integrals in form C - iS 

LOCAL VARIABLES 


Name 

Type 

Description 

n 

C 

Intermediate variables from Abromovitz and Stegun 

H 

C 

To rotate G in complex plane. 


FUNCTIONS 

1. To calculate a value for the Fresnel integrals 

SUBPROGARAMS CALLED 
None 

CALLING PROGRAMS 
NOPLFT 

ERRORS 

NONE 

ENTRY 

Calculate G from Abromovitz and Stegun equations p 302. 
Calculate H. 

Calculate E = C - iS following Abromovitz and Stegun. 
EXIT 



* ROUTINE - NQPLFT 

* 

* PURPOSE - To calculate the airfoil gust response function. 

* 

* AUTHOR - Roy X. Amiet 

* 

* INPUT 

* 

* ARGUMENT 

* 


* 

* 

Name 

Type 

Decsr iption 


X 

RD 

Chordwise distance 

* 

SG 

RD 

Far-field distance modified by Prandtl-Glauert 

* 



factor 

* 

RM 

RD 

Chordwise Mach number component 

* 

XX 

RD 

Chordwise wavenumber 

* 

ft 

YK 

RD 

Spanvise wavenumber 

t 

ft 

OUTPUT 



* 

ft 

ARGUMENT 



it 

* 

Name 

Type 

Description 

* 

GL2 

RD 

Square of affective lift, including 

* 

ft 



noncompactness 

* 

* 

LOCAL VARIABLES 


* 

* 

Name 

Type 

Description 

* 

B2 

RD 

Column vectors for which the cross product is 

* 



found. 

it 

UM 

RD 

Normalizes the output to a magnitude of 1. 

it 

RMI 

RD 

Similarity Mach number for skewed gusts 

it 

UMI 

RD 

Specifies use of low or high frequency solution 

it 

T2 

RD 

Intermediate dummy variable. 

t 

CM 

RD 

Magnitude of E. 

* 

E 

RD 

Complex representation of Fresnel integrals; 

* 



C - iS 

t 

FUNCTIONS 



t 

1. To calculate the 

airfoil response function for given values of 

it 

ft 

the vavenumbers 


it 

SUBPROGARAMS 

CALLED 


* 

ft 

NOPFNL 



t 

CALLING PROGRAMS 


* 

NOPPX 




* 


* ENTRY 

* Calculate Prandtl-Glauert factor 32. 

* UM is proportional to airfoil chord/acoustic wavelength. 

* RMI i3 Graham ' 3 similarity Mach number for skewed gusts. 
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UMI used in following line for branching test 
Case of (UMI -.75) HEG, ZERO, POS 
8EG: Use modified Sear3 function 
ZERO: Use modified Sears function 


POS: Use high reduced frequency solution 

Call Fresnel integrals 
Find absolute value 
Calculate square of effective lift 
End case. 

EXIT 
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* ROUTINE - MOPHI 


* PURPOSE - To taka tha cro33 product of column vectors M2 and 113 of 

* E and place the result in column 111 

* 


* 

t 

• 

* 

* 

* 

• 

* 

AUTHOR - Roy 
INPUT 
ARGUMENT 
Name 

K. Amiet 

Type Description 

• 

CL 

RD 

Chord/turbulence integral length scale 

* 

CVT 

RD 

Step size for summation over wavenumber 

♦ 

XK 

RD 

X wavenumber component 

* 

YK 

RD 

Y wavenumber component 

* 

ZK 

RD 

Z wavenumber component 

* 

ALP 

RD 

Angle alpha 

* 

G 

RD 

Gamma 

• 

DXDZ 

RD 

Deformation tensor; 9 components 

* 

* 

OUTPUT 



* 

ARGUMENT 



it 

it 

Name 

Type 

Description 

* 

SVV 

RD 

Pressure averaged around azimuth 

* 

NP 

I 

Counter for upward summation in NTITRB 

* 



subroutine 

• 

NM 

i. 

Counter for dovnvard summation in NTITRB 

* 



subroutine 

* 

N 

I 

Counter for rescaling of CVT in NTITRB 

* 



subroutine 

* 

* 

LOCAL VARIABLES 


* 

# 

Name 

Type 

Description 

• 

EX 

RD 

Constant in Karman spectrum 

• 

DK 

RD 

Wavevector; 3 components 

it 

Cl 

RD 

Cosine alpha 

it 

SI 

RD 

Sine alpha 

• 

C2 

RD 

Cosine gamma 

• 

S2 

RD 

Sine gamma 

it 

VN 

RD 

Normal to blade 

it 

ROT 

RD 

Expression in analysis 

# 

DXDZR1 

RD 

Intermediate deformation matrix 

it 

DXDZR 

RD 

Final deformation matrix after rotation 

it 

CVT1 

RD 

Dummy for CVT 

It 

ED 

RD 

Downstream unit vectors parallel to axes 

t 

EU 

RD 

Upstream unit vectors parallel to axes 

it 

DKM 

RD 

Wavevector amplitude downstream at rotor face 

it 

DKUK 

RD 

Ratio of downstream to upstream wavevector 

it 



amplitudes 
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UQDQ 

RD 

Ratio of upstream to downstream 
amplitudes 

velocity 

XKU0 

RD 

Upstream x value of vavevector; 

initial value 

YKU0 

RD 

Upstream y value of vavevector; 

Initial value 

XY0 

RD 

Expression in Karman spectrum; 

initial value 

DKM1 

RD 

Wavevector amplitude downstream 

at rotor face 

DRUK1 

RD 

Ratio of dovnstream to upstream vavevector 
amplitudes 

UQDQ1 

RD 

Ratio of upstream to dovnstream velocity 
amplitudes 

XKU 

RD 

Upstream x value of vavevector 


YKU 

RD 

Upstream y value of vavevector 


ZKU 

RD 

Upstream z value of vavevector 


XY 

RD 

Expression in Karman spectrum 


DH 

RD 

Denominator in Karman spectrum 


DELI 

RD 

Temporary variable for testing 

end of iteration 

AN 

RD 

Float N 


ZK1 

RD 

Value of x wavenumber 


ZK2 

RD 

Value of y wavenumber 



SUBPROGARAMS CALLED 
NOPKVC 

CALLING PROGRAM 
NOPPX 

ERRORS 

NOME 

ENTRY 

Calculate ke which equals 0.7468342/L equals 0.373417(c/L)/b; 
Define vavevector in airfoil fixed coordinate system; x along 
chord, y along span, and z along normal. 

Sine and cosine of angles alpha and gamma. 

Calculate normal to airfoil in 3ame coordinate system as for DR. 
Calculate rotation matrix. The first three components give the 
chord direction, the second three give the spanvise 
direction and the third three give the normal direction. 

DO 50 1=1,3 
DO 30 J=1 ,3 

Initialize DXDZR1 
DO 90 K=i,3 

Calculate DXDZR1 
END DO 
END DO 
END DO 
DO 70 1=1,3 
DO 60 J=1 ,3 

Initialize DXDZR 
DO 91 K=1 , 3 

Calculate DXDZR 
END DO 
END DO 
END DO 

Introduce dummy for CVT 



Gall deformation program NOPKVC 

Calculate denominator for upstream turbulence with kz - 0. 

This will be used for comparisonlater in program. 

Make CVT as large a3 possible without losing accuracy. Increase 
by factors of two and test size of DELI. 

Initialize sum variable, SVV. 

Begin summation with minimum but still positive value of ZK after 
subtracting integer number of CVT values. 

Initialize counters HP and HM for number of steps up and down in kz. 
Sum over increasing values of ZK. 

Increment upward summation counter 
Set third wavevector component 
Call deformation program NOPKVC 

Calculate upstream values of three wavevector components 

Calculate XY, an intermediate variable 

Calculate DH, the denominator of the Karman spectrum 

Add to summation, SVV 

Check for excessive summations 

Increment z wavevector component 

Check to 3es how large the denominator ha 3 become. 

End upward summation 

Sum over decreasing values of ZK. 

Increment downward summation counter 
Decrement z wavevector component 
Set third wavevector component 
Call deformation program NOPKVC 

Calculate upstream values of three wavevector components 

Calculate XY, an intermediate variable 

Calculate DN, the denominator of the Karman spectrum 

Add to summation, SVV 

Check for excessive summations 

Check to 3ee how large the denominator ha3 become. 

End downward summation 
Multiply SVV by remaining factors 
EXIT 



RppendiH R: Computer Program for Noise Produced by 
Rotor-Turbulence Interaction 


A computer program implementing the algorithm discussed in the Theory Manual for Noise 
Calculation is given in the Computer Listings to follow. Computer Listing 1 is the main executive 
portion of the program; it prompts for inputs and prints the output. It also performs the integration 
over the rotor span. 

The inputs to the program, including the symbols used in the program, are: 

CL c/L = chord/turbulence length scale 

BN Blade number 

C0C c 0 /c = sound speed/chord 

RC Blade radius/chord 

R ic r/c = far-field distance/chord 

RPS Revolutions/sec 

ZM M z = axial Mach number 

FM Flow Mach number in rotor plane 



TH 0 = observer angle wrt z axis; observer in x,y plane 

SUMS Number of terms in azimuthal summation; N in equation (B3) 

PSI i|i = angle of M f wrt z axis; M f is positive inward as in figure 1 

DH Frequency step size in fractions of a harmonic 

AL1 Initial harmonic number 

L2 Number of frequency steps 

The numerical constant 203.95 occurring in the expression for PSD is obtained by introducing 
equation (82) into equation (63) and multiplying by the blade number. This gives 


S pp( ^o, x) 



/to Zsppb \ 2 

V c 0 a 3 2 * 


tt U b d |£| 2 b 2 F 67 


where 


F s 


2rr n 

b 2 u 2 Z “ 
n = 


^ ww (K x ,Ky,K 2 ^) 


(Al) 


Atmospheric pressure P 0 ■ 10 6 /0.987 dynes/cm 2 and reference pressure P r = 0.0002 dynes/cm 2 . 
Spp must be multiplied by 4tt to convert from a two sided - «<(o<ootoa one sided 0 < u < ® 

spectrum and to convert to a unit Hertz rather than a unit radian frequency of measurement. Using the 
relations 

K x = cob/U b c 0 2 = 7P 0 /Po (A2) 

this can be written (where 7 = 2tto/N) 
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Qpp(fo)/Pr J = (BKx) 2 lil^F} (A3) 


Converting to dB, 


10 log 


[iee)] — 


10 


(A4) 


The outputs of the program are: 

H Frequency of sound expressed in multiple of blade passage harmonic 

FI Frequency of sound in Hz 

PSD Spectrum level in dB (per unit Hz) 

NP 1,NM 1,N Counters for certain of the summations for use in diagnosis of output 


The programs in the following appendices are standard Fortran code with added line numbers for 
reference in the documentation. This allows more detailed comments to be used without deleter iously 
affecting the ability to easily read the code. 

The program listing here is for an homogenous turbulence in the rotor plane. This can be easily 
extended to the nonisotropic case by inputting the deformation tensor at each point of the rotor plane 
rather than just once. The turbulence spectrum at each point of the rotor plane is then calculated from 
a rapid contraction of an isotropic turbulence using this deformation tensor. The changes to the main 
program and one subroutine for this generalization to the nonhomogeneous case are listed below. 


Changes to main program: 
Line * Cfldfi 


12 

1 


RADIUS/C 

:, FAR-FIELD/C, RPS, IH ' /) 


24 

IF 

(IH 

.EQ. 0) 

READU,*) ( ( DXDZ ( I , J) , 1 = 

1,3), J= 

25 

IF 

(IH 

.EQ. ,0) 

WRITE(* ,200 ) ( (DXDZ ( I , J) , 

J=1 ,3) , 

26 

IF 

(IH 

.EQ. 0) 

WRITE(2,200 ) ((DXDZ(I.J), 

J=1 ,3) , 

28 

IF 

(IH 

.EQ. 0) 

CLOSE ( 1 ) 


59 

1 



P2,NP1 ,NM1 ,N, IH) 



.3) 

1=1,3) 

1=1,3) 


In the subroutine SPX, change the first line of the subroutine to accept the new variable IH and add the 
following line between lines 46 and 47: 


1 SUBROUTINE PX ( F 1 , ZM , Cl , BN , C0C , UU , RC , R 1C , RPS , FM , TH , J , PS I , DXDZ , 

46a IF (IH .EQ. 1) READU,*) (CDXDZ(I.J), J=l,3), 1 = 1,3) 


In the above code the parameter IH is an input designating whether a homogeneous or inhomogeneous 
calculation should be performed. If IH = 0 is input, homogeneous turbulence is assumed. If IH = 1 is 
input, inhomogeneous turbulence is assumed, and the above read and write statements are deferred to 
the subroutine SPX. 
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1 

2 

3 

4 


5 

500 

6 

1 

7 


8 


9 


10 

30 

11 

100 

12 

1 

13 


14 


15 


16 


17 

300 

18 

1 

19 

1 

20 

1 

21 

10 

22 


23 


24 


25 


26 


27 

200 

28 


29 

20 

30 

400 

31 


32 


33 


34 

900 

35 

1 

36 

70 

37 

700 

38 

1 

39 

1 

40 

1 

41 


42 


43 


44 


45 


46 


47 


48 

800 

49 


50 



PXFINI 12/3/84 R. K. Amiet Integrates over, span. 

DIMENSION 0XDZC3.3) 

CHARACTER* 10 FNAME 
CHARACTER* 10 FNAME 1 
WRITE (* ,500 ) 

F0RMAT( ' Turbulence ingestion program with integration over' 

span'/' and nonisotropic turbulence. R. K. Amiet 12/84'/) 
WRITEC*, ' (A\) ' ) ' OUTPUT FILE NAME 
READ(* , ' (A)' ) FNAME 1 
OPEN (2 . FILE=FNAME 1 , STATUS= ' NEW ' ) 

WRITEC*, 100) 

FORMAT (//' INPUT C/L, BLADE NO., C0/C, ' 

' RADIUS/C, FAR-FIELD/C , RPS'/) 

REAOC*,*) CL, BN, C0C , RC , R1C, RPS 

TM1 = 2. *3. 141593*RPS*RC/C0C 

WRITEC*, 300) TM1, CL, BN, C0C, RC, R1C, RPS 

WRITEC2.300) TM1, CL, BN, C0C, RC, R1C, RPS 

FORMAT C//‘ M TIP = ',F4.3,' CHORO/TURB SCALE =',F7.4/ 

' BLADE NO. = ' , F4 . 0 , ' SOUND SPEED/CHORD =',F7.1/ 

’ RADIUS/CHORD =',F5.1,' FAR-FIELD/CHORD =', 

F6. 1/' REV/SEC = ' ,F5. 1/) 

WRITEC*, ' CA\) ' ) ' INPUT FILE NAME -' 

REAOC*; 'CAD FNAME 

OPENC 1 ,FILE=FNAME ,STATUS= ' OLD ' ) 

REAOC 1 ,*) ( (DXDZC I , J ) , 1 = 1,3), J=l,3) 

WRITEC*. 200) (CDXDZCI.J), J=l,3), 1=1,3) 

WRITEC2.200 ) (CDXDZCI.J), J=l,3), 1=1,3) 

FORMAT ( 1X.3F5 .2 ) 

CLOSE ( 1 ) 

WRITEC*, 400) 

FORMATC/' INPUT MZ, FM, UU, THETA, SUMS, PHI'/) 

REAOC*,*) ZM, FM, UU, TH, J, PH 

WRITEC*, 900) ZM, UU, TH. FM, J, PH 

WRITEC2.900) ZM, UU, TH, FM, J, PH 

FORMATC MZ =' ,F4. 3,' UU=',F6.3,' THETA =' ,F5 . 0/ 

MF =' ,F4.3, ' SUMS = ' , 14, ' PSI=',F5.1/) 

WRITEC*, 700) * 

FORMATC/' INPUT'/ 

BEGINNING HARMONIC: If < 0, go to matrix input.'/ 

FRACTIONAL HARMONIC SPACING: If < 0, go to previous input.'/ 

NUMBER STEPS: If < 0, go to first input.'/) 

READC*,*) AL1 , OH, L2 
IF CALI .LT. 0.) GO TO 10 
IF (OH .LT. 0.) GO TO 20 
IF (L2 .LT. 0) GO TO 30 
DF = RPS*BN*DH 
WRITEC*, 800) 

WRITEC2 ,800 ) 

FORMAT ( /2X , ' HRMNC ' , 5X , ' FREQ ' , 6X , ' PSO ' , 4X , ' SM P',1X,'SM M'/) 

FI = CAL 1/OH- 1. )*DF 
IF (L2 .EQ. 0) L2= 1 


Computer Listing 1: Executive program for turbulence Ingestion. 
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51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 


DO 60 K= 1 , L2 
FI = FI + DF 
RCP = RC*1.05 
SUM * 0. 

DO 22 12=0,19 
RCP = RCP - . 05*RC 
TM = TMPRCP/RC 

CALL PX (F 1 , ZM ,CL , BN , C0C ,UU ,RCP ,R 1C , RPS ,FM , TH , J,PH ,DXDZ , 
1 P2 , NP 1 , NM 1 , N ) 

IF (12 .EQ. 0) P2 = P2/2 . 

SUM = SUM + P2 
22 CONTINUE 

PSD = 10 . *ALOG 10 (SUM*RC* . 05 ) 

H = F1/(RPS*BN) 

WRITEC* ,600 ) H, FI, PSD, NP1, NM 1 , N 
WRITEC2 ,600 J H, FI. PSD, NP1, NM 1 , N 
600 FORMAT ( IX, F6. 1, IX, F8. 1, IX, F9. 2, IX, 15, 216) 

60 CONTINUE 
GO TO 70 
END 


Line * Comment 

7-9 Prompt for a filename into which to write the calculated results 
10-13 Prompts for inputs 
14 Calculation of tip Mach number 

15-20 Prints the input quantities for checking 

21-23 Prompt for input containing the deformation tensor values 
24-28 Read and write the values for the deformation tensor 
29-35 Prompt for more input parameters; writes out values to file 
36-44 Prompt for final input values; if any are negative the program returns to one of the 
previous input prompts 
45 DF = frequency step size 

46-48 Prints headings for the eventual printout 

49 Initialize the frequency variable 

50 Input of zero for L2 gives one iteration with this line 

51-69 Stepping through frequency range 

52 Increment frequency variable 

53-54 Initialize radius and summation variables 

55-63 Integration over radius 

56 RCP = local radius 

57 TM = M t local tip Mach number 

58-59 Call routine for summation over 7 

60 Divide the first term in the summation by 2; the last term also should halved, but since it 
approaches zero, it needn't be explicitly set to zero. 

61 Perform radial summation 

63 Introduce remaining factors inequation (69) 

64 Calculate harmonic number 

65-67 Write outputs 


Computer Listing 1: Executive program for turbulence Ingestion. 
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Appendix B: Subroutine Calculating the Main Problem Parameters 

and Integrating ouer Azimuth 

This subroutine is called by the main turbulence ingestion program in Appendix A. It calculates 
geometry and integrates over azimuthal angle. The inputs to the subroutine are: 


FI Observer frequency 

ZM Axial Mach number; see figure 1 

CL Chord/turbulence integral length scale 

BN Blade number 

COC Sound speed/chord 

UU RMS turbulence intensity/axial velocity 

RC Local blade radius/chord 

R1C Far-field distance/chord 

RPS Revolutions per second 

FM Flight Mach number in rotor plane; see figure 1 

TH Polar angle of observer; see figure 1 

J Number of azimuthal integration points 

PSI Angle of flight Mach number with the y axis; see figure 1 

DXDZ Deformation tensor 


The outputs are: 

P2 Sound level from a radial blade segment 

NP1 ,NM1 ,N Counters used in the summations; used for diagnostic purposes 

The program calls the subroutines TRBNI and LEFF, the turbulence and the airfoil response 
calculations. The subroutine performs the azimuthal integration over © in Eq. (63). 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 


• Subroutine SPX 

• Program written by R. K. Amiet. 

SUBROUTINE PX(F1 , ZM,CL,BN,C0C,UU,RC,R1C,RPS,FM,TH, J,PSI ,DXDZ, 
1 P2,NP1,NM1,N) 

DIMENSION 0XDZ(3,3) 

PI = 3.14159 

TM = 2 .*PI*RPS*RC/C0C 

AJ = J 

DEL = 360. /AJ 
C = COSCTH/57 .2958) 

S = SIN(TH/57.2958) 

C4 = C0S(PSI/57 .2958) 

54 = SIN (PS 1/57. 2958) 

QM = FM*S*S4+ZM*C 
BFZ2 = l.-FM**2-ZM**2 
SQ1 = SORT (QM**2+BFZ2 ) 

RER = (SQ1+QM)/BFZ2 

NP 1 = 0 
NM1 = 0 
SUM = 0. 

G = 90. - DEL 
DO 50 1=1, J 
G = G + DEL 
C2 = COS (G/57. 2958) 

S2 = SIN(G/57 .2958) 

C5 = C0S( (G+PSI )/57 .2958) 

55 = SIMC (G+PSI )/57. 2958) 

ALP = 57 . 2958*ATAN2 ( ZM , TM+FM*C5 ) 

Cl = C0S(ALP/57 .2958) 

SI = SIN (ALP/57. 2958) 

RM = SQRT (ZM**2 + (TM + FM*C5)**2) 

C3 = C*S1 - S*C 1*S2 
X = R1C*(TM*RER*C1 - C3) 

Y = R1C*(S*C2 + FM*RER*S5) 

Z = R1C*(C*C1 + S*S1*S2 + TM*RER*S1 ) 

SG = SQRT (X**2,+ ( 1 .-RM**2)*(Y**2+Z**2) ) 

FP = F1*(1.+TM*(S2*S-FM*RER*C5)/SQ1) 

XK = PI*FP/(C0C*RM) 

YK = PI*FP*Y/(SG*C0C) 

T = 2.*PI*RC/(TM*C0C*BN) 

T 1 = T*TM*S1*C1/ZM 

• Next line added to make peaks occur at bpf when ZM = 0 . 

T 1 = T/(l.+ MF*C5/TM) 

XX = TM* (T-T 1 ) 

YY = T 1*FM*S5 
ZZ = T*TM*C0C*S1 

T2 = T 1 + XX*(RM-X/SG)/( 1 .- RM**2) + YY*Y/SG 

CVT = PI/ZZ 

ZK0 = PI*T2*F 1 /ZZ 

CALL TRBN I ( CL , CVT , XK , YK , ZK0 , ALP , G , DXDZ , SVV , NP , NM , N ) 


Computer Listing 2: Subroutine for calculating principle parameters in problem 


24 


48 CALL LFT (X,SG,RM,XK,YK,GL2) 

49 NP1 = NP1 + NP 

50 NM1 = NM1 + NM 

51 AD = GL2*SW*RM*(RM*Z*FP*XK/(F1*SG**2) )**2 

52 SUM = SUM + AD 

53 50 CONTINUE 

54 P2 = 2.48E20*SUM*BN*(UU*ZM)**2/(C0C*AJ1 

55 RETURN 

56 END 

Line # Comment 

1-2 Inputs 

5 TM = M t = local tip Mach number 

6 Float J 

7 DEL = azimuthal step size 

8-11 C, S and C4, S4 = cosine and sine of 6 and i|r respectively 

12 QM = M $ cos©;Eq. (43) 

13 BFZ2 = 1 - M $ 2 ;Eq. (43) 

14 SQ1 = ( 1 - M $ 2 sin 2 ©) 1/2 

15 RER = r e /r;Eq. (42) 

16-19 Initialize variables 

20-53 Integration over azimuthal angle y 

21 Increment r 

22-25 C2, S2 and C5, S5 = cosine and sine of y and y + # respectively 

26 ALP =a;Eq. (34) 

27-28 Cl, SI = cosine and sine of a 

29 RM = M b ;Eq. (36) 

30 C3 = cos $; Eq. (49) 

31-33 X,Y,Z = x 3 , y 3 ,z 3 ;Eq.(48) 

34 SG = a 

35 FP = frequency on blade; Eq. (55) 

36-37 XK, YK = specific wavenumbers K x , Kyj Eqs. (9), ( 14) 

38 T = time between blade passes 

39-40 T 1 = time between eddy intersections; Eq. (62) for line 39. Equation (32) is not 

appropriate for eddies stretched in axial direction; for this case the diagonal distance in 
figure 3 is not appropriate; given here is the result for an eddy of infinite axial length 
corresponding to the limit M z -* 0; Eq. (62). 

41-43 XX, YY, ZZ = X, Y, Z in Eqs. (62), (64) 

44 CVT = 2tt/Z = step size in Eq. (70) 

45 ZK0 = cooTVZ = initial radian frequency in Eq. (70) 

46 T2 = T 2 ; Eq. (65) 

47 Call turbulence summation subroutine 

48 Call airfoil response subroutine 
49-50 Increment counters 

51 Contribution to spectrum from azimuthal integration; Eqs. (50) and (69) 

52 Summation over azimuthal spectrum 

54 Multiplication by remaining factors in Eq. (69) 


Computer Listing 2: Subroutine for calculating principle parameters in problem 
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Appendix C: Subroutine for Summing ouer the Turbulence Spectrum for 

Blade to Blade Correlation 

This subroutine calculates the summation in Eq. (69) over the turbulence spectrum. The output 
SVV of this program is 


SVV = 


2tt 

b 2 u 2 Z 


Z * WW (X,Ky,K ( z n) ) 


(El) 


The inputs to the program and the equivalent symbols used in the report are 

CL = c/L CVT = 2tt b/Z 

XK = K x YK = K y 

ZK0 = co 0 T 2 /Z DXDZ(I,J) = (x^ = precontraction coordinate; = post 

contraction coordinate) 

Other parameters in the program are: EK = k e b where k e is given in Eq. (32). CVT is the step size in 
the summation. Since the spectrum is written in a form with the wavenumbers k x and k y normalized 
by k e , CVT must also be so normalized. Note that this changes CVT in the main program, but this is of 
no concern since CVT is used only in this subroutine. 

The summation over n is over the range ±oo. in numerically carrying out this summation it is 
easiest to begin where the individual terms are maximum and sum to either side of this point until the 
terms become negligible. Thus, the integer N is the optimum value of n at which to begin the 
summation. The optimum value would be n = 0 except for the offset <o 0 T 2 /Z occurring in K z ^ n ^ (see 

Eq. (69)). 

If the step size CVT is too small, very many terms in the summation would be required. For this 
case, the summation can be replaced by an integral. For the case of an isotropic turbulence the 
integral can be carried out in closed form as in Eq. ( 3 1 ). 

The parameters ZK 1 and ZK2 represent K z ^ n V k e . ZK 1 is used for increasing values of n and 
ZK2 for decreasing values. 

The integers NP and NM count the number of summation steps in the upward and downward 
directions respectively. If either goes over 100 the summation is terminated, although this case has 
never been encountered. If the output specifies that NP or NM reached the value of 100, further 
checking should be done to ensure that a sufficient number of steps were taken. Generally the 
summation is terminated when the term DN (denoting the denominator in Eq. (32) ) becomes greater 
than 200 times the initial value of the denominator. 

Finally, the coefficient 0.1 1561 = 55 r(5/6)/[9/rT T( 1/3) 4tt] . The factor 4rr comes from 
Eq. (30) and the remaining factors come from the expression for I. 
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EDDYNI.FTN; Summation over turbulence spectrum. 
Written by R. K. Amiet; final version 1985 

3N I ( CL , CVT , XK , YK , ZK , ALP , G , DXDZ . SVV , NP , Nl 
VM(3) , DKC3), DXDZ (3,3) , 

, ROT (3,3), 0X0ZR1 (3,3) 


1 


SUBROUTINE TRBNI (CL ,CVT 

2 


DIMENSION EU (3,3) , ED(3,: 

3 

1 

DXDZR (3, 

4 


EK = 373417*CL 

5 


DK( 1 ) = XK 

6 


DK(2) = YK 

7 


DK(3) = 0. 

8 


Cl = COS ( ALP/57 . 2958 ) 

9 


SI = SIN(ALP/57.2958) 

10 


C2 = C0S(G/57 .2958) 

11 


S2 = SIN(G/57 .2958) 

12 


VN( 1 ) = 0. 

13 


VN(2) = 0. 

14 


VN(3) = 1. 

15 


ROT (1,1) = C1*S2 

16 


ROT (2,1) = C2 

17 


ROT (3,1) = S1*S2 

18 


ROT (1,2) = -C1*C2 

19 


ROT (2,2) = S2 

20 


ROTO, 2) = -S1*C2 

21 


ROT (1,3) = -SI 

22 


ROT (2,3) = 0. 

23 


ROTO, 3) = Cl 

24 


DO 50 1=1,3 

25 


DO 30 J=1 ,3 

26 


DXDZR 1 ( I , J ) = 0. 

27 


DO 90 K=1 ,3 

28 


DXDZR 1 ( I , J) = DXDZ(I,K)*I 

29 

90 

CONTINUE 

30 

30 

CONTINUE 

31 

50 

CONTINUE 

32 


DO 70 1=1,3 

33 


DO 60 J=1 ,3 

34 


DXDZR ( I , J) = 0. 

35 


DO 91 K=1 ,3 

36 


DXDZR ( I , J) = DXDZRKK, J) 

37 

91 

CONTINUE 

38 

60 

CONTINUE 

39 

70 

CONTINUE 

40 


CVT 1 = CVT 

41 


CALL KVEC(DK, DXDZR, VN, ED 

42 


XKU0 = EU( 1 ,3)*DKM/DKUK 

43 


YKU0 = EU(2,3)*DKM/DKUK 

44 


XY0 = EK**2 + XKU0**2 + 

45 


D0 = XY0**2. 83333 

46 


N = 0 

47 

10 

N = N + 1 

48 


CVT = CVT 1 


OXDZRKI , J) 


DXOZRC I , J) 


YKU0**2 


Computer Listing 3: Subroutine for calculating summation over turbulence spectrum. 
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49 


IF (N 

.GT. 1) CVT 1 = 2 .*CVT 1 

50 


DK(3) 

= CVT 1 

51 


CALL 

KVEC ( DK , DXDZR . VM . ED , EU , DKM 1 . DKUK 

52 


XKU = 

EU ( 1 , 3 ) *OKM 1 /OKUK 1 

53 


YKU = 

EU (2,3) *OKM 1 /DKUK 1 

54 


ZKU = 

EU ( 3 , 3 ) *OKM 1 /DKUK 1 

55 


XY = 

EK**2 + XKU**2 + YKU**2 

56 


DM = 

(XY + ZKU**2)**2. 83333 

57 


DELI 

= ABS( 1 .- OM/O0 ) 

58 


IF (DELI .LT. .01) GO TO 10 

59 


SVV = 

0. 

60 


N = ZK/CVT 

61 


AM = 

M 

62 


ZK1 = 

ZK - AM*CVT 

63 


ZK2 = 

ZK1 

64 


MP = 

0 

65 


MM = 

0 

66 

20 

MP = 

MP + 1 

67 


DK(3) 

= ZK1 

68 


CALL 

KVEC ( DK . OXOZR . VM , ED , EU , OKM . DKUK , 1 

69 


XKU = 

EU( 1 ,3)*DKM/DKUK 

70 


YKU = 

EU ( 2 , 3 ) *DKM/OKUK 

71 


ZKU = 

EU ( 3 , 3 ) *DKM/DKUK 

72 


XY = 

EK**2 + XKU**2 + YKU**2 

73 


DM * 

(XY + ZKU**2) **2.83333 

74 


SVV = 

SVV + (XY - EK**2 ) / ( UQDQ**2*DM ) 

75 


IF (MP .GT. 100) GO TO 40 

76 


ZK1 = 

ZK1 + CVT 

77 


IF (DM .LT. 00*200.) GO TO 20 

78 

40 

MM = 

MM + 1 

79 


ZK2 = 

ZK2 - CVT 

80 


DK(3) 

= ZK2 

81 


CALL 

KVEC ( DK , DXOZR . VM . ED , EU . DKM , OKUK . 

82 


XKU = 

EU ( 1 , 3 ) *DKM/OKUK 

83 


YKU = 

EU ( 2 , 3 ) *DKM /OKUK 

84 


ZKU = 

EU ( 3 , 3 JTDKM/DKUK 

85 


XY = 

EK**2 ♦ XKU**2 + YKU**2 

86 


OM = 

(XY +ZKU**2) **2.83333 

87 


SVV = 

: SVV * (XY - EK**2)/(UQ0Q**2*0M) 

88 


IF (MM .GT. 100) GO TO 80 

89 


IF (OM .LT. 00*200.) GO TO 40 

90 

80 

SVV = 

> . 1 156 1*SVV*CVT*EK** . 666667 

91 


RETURN 

92 


EMD 
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Line* Comment 

1 Inputs 

4 k e = 0.7468342/L = 0.373417(c/L)/b; 0.7468342 = r(5/6)/r T /r( 1/3). 

5-7 Define wavevector k in coordinate system fixed to airfoil; x along chord, y along span, 

and z along normal. k z will take different values. 

8-11 Sine and cosine of angles a and y. 

12-14 Normal to airfoil in same coordinate system as in lines 5-7. 

15-23 Rotation matrix; dx t /dxz with the definitions in Eq.( 47 ). The first three components 
give the chord direction, the second three give the spanwise direction and the third three 
give the normal direction. 

2 1 VN( 1 ) written in the unrotated x t system. 

22 VN(2) written in the unrotated & system. 

23 VN(3) written in the unrotated Xj system. 

24-39 Calculate deformation tensor in coordinate system defined in lines 5-7. First 
premultiply, then postmultiply, DXD2 by ROT. 

40-45 Denominator for upstream turbulence with k z = 0. Used for comparison in lines 77 and 
89. 

46-58 If CVT is too small, will take too many summations. Make CVT larger until start getting 
significant variation in either the ratio of wavevectors between downstream and 
upstream. 

59 Initialize sum variable. 

60-63 Begin with minimum but still positive value after subtracting integer number of CVT 
values. 

64-65 Initialize counters for number of steps up and down in k r 

66-77 Sum over increasing values of k z . Lines 69-71 calculate upstream values of 

wavevector. These expressions are placed into isotropic turbulence in 72-83. Line 75: 
check on number of iterations. Line 76: increment k z . Line 77: check to see how large 

the denominator has become. 

78-89 Same comments as lines 53-64 except summation is in downward direction. 

90 Factors multiply sum. Note that CVT is included in the product so that it is o'k to 

increase CVT in lines 46-58 making the summation step size bigger. 

0.11561 = 55r(5/6)/[36TT 8/2 r( 1/3)] 
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Appendix 0: Program for Calculating Airfoil Response Function, £, 

for Leading Edge Noise 


Computer Listing 4 gives the subroutine for calculating the effective lift |£| for noise produced 
by airfoil- turbulence interaction. The inputs to the subroutine and the equivalent notation used in 
this report are: 


Program 

Text 

X 

*3 

SG 

a 

RM 

M b ; equation (36) 

XK 

k x ; x wavenumber component 

YK 

ky ; y wavenumber component 


The output 6L2 is the effective lift |£| 2 and is calculated in a low frequency and a high frequency 
regime. The low frequency result is given by equation (27) and the high frequency result by Eqs. 


(28). 



* 


Lift function subroutine 

* 


Written by R. K. Amiet; find version 1978 

1 


SUBROUT I ME LFT (X,SG,RM,XK,YK,GL2) 

2 


COMPLEX E 

3 


62 = 1 . -RM**2 

4 


UM = RM*XK/B2 

5 


RMI = SORT (RM**2-B2*(YK/XK)**2) 

6 


UMI = UM*RMI/RM 

7 


IF (UMI- .75) 10, 10, 20 

8 

10 

0L2 = 1 . / ( B2* ( 1 . / ( 1 . +2 . 4*XK/B2 ) +6 . 283 19*XK/B2 ) ) 

9 


GO TO 30 

10 

20 

T 1 = ABS(UM*(X/SG-RMI/RM) ) 

11 


CALL FRMLCSQRT ( 1 .27324*T 1 ) ,E) 

12 


CM = CABS(E) 

13 


GL2 = . 20264*CM # *2/(XK*T 1*( 1 . +RMI ) ) 

14 

30 

RETURM 

15 


END 

Line * 


CflmmfiDis 

1 


Inputs and output 

3 


B2 * 8 b 2 

4 


UM = u 

5 


RMI = M.; equation (24) 

6 


UMI = u.; equation (24) 

7 


Check whether to use low or high frequency solution. 

8 


Low frequency solution; equation (27) 

10-13 


High frequency solution; equation (28a) 

10 


T 1 = ©i; equation (29) 

13 


Equation (28a); .20264 = 2/tt 2 
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Appendix E: Program for Calculating the Effect of Contraction on 

the Turbulence Spectrum 

Computer Listing 5 gives the subroutine for calculating the effect of a contraction of a 
turbulence spectrum. This subroutine has as input a 3x3 deformation matrix representing a flow 
contraction. The program then accepts a flow direction V n and a wavevector k both for the downstream 

flow, and calculates the upstream wavevector and flow vector. The inputs to the subroutine and the 
equivalent notation used in this report are: 

Program, Iaai 

DK k d = downstream wavevector components 

DXDZ dXj/5£j = deformation tensor; equation (71) 

VN Upwash velocity component normal to the airfoil 

ED, EU Tensors representing the three vector components of the upstream and downstream 
coordinate systems. Each such vector is represented by a column of the matrix. 


The outputs are: 


DKM 


Wavevector amplitude downstream at the rotor face 

DKUK 


Ratio of downstream to upstream wavevector amplitudes 

UQDO 


Ratio of upstream to downstream velocity amplitudes 

* 


Turbulence contraction calculation. 

* 


Written by R. K. Amiet 9/7/83 

1 


SUBROUT I ME K VEC ( DK , DXDZ , VM , ED , EU , DKM , DKUK , UODQ ) 

2 


DIMEMSIOM EDO, 3), EUO,3), VMC3) , F3(3), DK ( 3 ) , OXOZO.3) 

3 


DKM = SORT (DK( 1 )**2 + DK(2)**2 + 0K(3)**2) 

4 


DO 10 1-1,3 

5 


ED (1,3) = OKCD/DKM 

6 


ED (1,2) = VM( I ) 

7 

10 

COMTIMUE 

8 


CALL CROSS (ED, 1 ) 

9 


CALL CROSS (ED, 2) 

10 


DO 20 1=1,3 

11 


DO 30 J=1 ,3 

12 


EU ( I , J ) = 0. 

13 

30 

COMTIMUE 

14 

20 

COMTIMUE 

15 


DO 40 I-l,$ 

16 


DO 50 J=1 ,3 

17 


DO 60 K=1 ,3 

18 


EU ( I , J ) = EU( I , J) + DXDZ(K, I )*ED(K, J) 

19 

60 

COMTIMUE 

20 

50 

COMTIMUE 

21 

40 

COMTIMUE 

22 


F3( 1 ) = EU (1,3) 

23 


F3(2) = EU(2,3) 

24 


F3(3) = EU(3,3) 

25 


AF1 = SQRT ( EU (1,1 )**2 + EU(2,1)**2 + EU(3,1)**2) 


Computer Listing 5: Subroutine for calculating effect of contraction on turbulence 
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26 

EU (1,1) = EUU/n/AFl 

27 

EU (2,1) = EU(2, D/AF1 

28 

EU (3,1) = EUC3,1)/AF1 

29 

CALL CROSS (EU, 3) 

30 

CALL CROSS (EU, 2) 

31 

DKUK = EU ( 1 ,3)*F3( 1 ) * 

32 

UQOQ = AF 1*DKUK 

33 

RETURN 

34 

END 

Jne * 

Comments 

3 

Calculates the magnitude of the 

4-7 

Reads V n into the 2'nd column o 


EU(2,3)*F3(2) + EU(3,3)*F3(3) 


Calculates the direction of the vorticity determined by the specified wavevector k and 
velocity V n directions. 


9 Calculates the direction of the velocity produced by the wavevector. V n is a component of 

this velocity. 

10-14 Initialize matrix e u . 

15-21 Calculates the values of e u from the downstream values e d by using equation(7 1 ) on each of 
the 3 column vectors making up e d . This is not the final form for e u . The vector e t u 
representing the vorticity direction will keep its direction but be normalized to a unit 
vector. e 3 u will be found by taking the cross product of vectors 1 and 2, and e 2 u is 
then determined by taking the cross product of e 3 u and e! U . 

22-24 F represents a vector found by taking a unit vector along k d and transforming according to 
equation(B 1 lb). This is used to determine k d /k u using equation(72). 

25 The magnitude of the transformed value of the vorticity, needed in equation(73). 

26-28 Unit vector in the direction of the upstream vorticity. 

29 Determines the unit vector e 3 u from the cross product of ej u and_f 2 u . 

30 Determines the unit vector e 2 u from the cross product of e 3 u and ei u . 

31 Determines k d /k u using equation( 72). 
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Appendix F: Subroutine for Calculating the Cross Product of Tuio Uectors 

Computer Listing 6 gives the subroutine for calculating the cross product of two vectors. The 
subroutine takes a 3x3 matrix E(i,j) composed of 3 column vectors ej and replaces the Nl'th vector 

by the normalized cross product of the other two in proper cyclic order .The inputs to the subroutine 
and the equivalent notation used in this report are: 

Program Iasi 

E(i,j) Tensor with columns composed of vectors; also represents the output 

N 1 Column not involved in the cross product 


* Subroutine for calculating cross product of two vectors 

* Written by R. K. Amiet 9/1/83 

1 SUBROUTINE CROSS(E.Nl) 

2 DIMENSION E(3,3) 

3 N2 = N1 + 1 

4 N3 = N1 + 2 

5 IF (N2 .GT. 3) N2 = N2 - 3 

6 IF (N3 .GT. 3) N3 = N3 - 3 

? E( 1 ,N1 ) = E(2,N2)*E(3,N3) - E(3,N2)*E(2,N3) 

8 E(2,N1 ) = E(3,N2)*E( 1 ,N3) - E( 1 ,N2)*E(3,N3) 

9 E(3,N1 ) = E( 1 ,N2)*E(2,N3) - E(2,N2)*E( 1 ,N3) 

10 EN1M = SORT (E( 1 ,N1 )**2 E(2,N1)**2 + E(3,N1)**2) 

11 E( 1 ,N1) = E( 1 ,N1 J/EN1M 

12 EC2.N1) = E(2,N1 1/EN1M 

13 E(3,N1) = E(3,N1)/EN1M 

14 RETURN 

15 END 

Lipe.* Comm, eat 

2-6 Defines the cyclic sequence 1,2,3 or 2,3,1 or 3,1,2 where N1 is the first number of 
the sequence. 

7-9 Replaces the Nl'th vector by the cross product of the other two. 

10-13 Normalizes the Nl'th vector to have a magnitude of 1. 
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Appendix G: Subroutine for Calculating the Fresnel Integral 


Computer Listing 7 gives the subroutine for calculating the Fresnel integrals. The program uses 
the algorithm given by equations 7.3.9, 7.3.32, 7.3.33 of Abromowitz and Stegun, Handbook of 
Mathematical Functions, Dover Publications, New York, 1968. 


* 


1 

2 

3 

4 

5 

6 

7 

8 


FRNL.FOR Subroutine for calculation of Fresnel integrals 
Written by R. K. Amiet 1976 
SUBROUTINE FRNL(X,E) 

COMPLEX E, G, H 

G = CMPLX ( ( 1 . + . 926*X ) / (2 . + 1.792*X + 3.104*X**2), 

1./C2.+ 4 . 142*X + 3 . 492*X**2 + 6.67*X**3)) 

H = CMPLX(SIN(1.5708*X**2),COS(1.5708*X**2» 

E = G*H + CMPLX ( .5, -.5) 

RETURN 

END 
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Appendix H: Sample Calculation 

These two test cases use the same inputs, listed below. The only difference between the cases is 
the deformation tensor. For the first case the tensor represents no deformation. The second 
represents a deformation by a factor of 4 in the axial direction. 

Rev/sec=10. Chord/Turb scale = .01 Machno.M f =0. 

Blade * = 2. Sound speed/Chord = 1000. Far field/chord = 100. 

Radius/chord = 10. rms turb vel/U = .0 10 Theta0 = 0. 

Axial Mach no. = M z = 0.1 

Case 1 

Deformation tensor = 

1.00 .00 .00 
.00 1.00 .00 
.00 .00 1.00 


HRMNC 

FREQ 

PSD 

HRMNC 

FREQ 

PSD 

.5 

10. 

.0 

31. 

.94 

5. 

.5 

110 

.0 

45.62 

1.0 

20. 

.0 

62, 

,82 

6. 

,0 

120 

.0 

45.58 

1 .5 

30. 

.0 

44. 

,56 

6. 

.5 

130 

.0 

44.66 

2.0 

40. 

.0 

56. 

.06 

7, 

.0 

140 

.0 

44.39 

2.5 

50. 

.0 

47, 

. 12 

7. 

.5 

150 

.0 

43.73 

3.0 

60. 

.0 

51. 

.93 

8, 

.0 

160 

.0 

43.36 

3.5 

70 

.0 

47 

.24 

8 

.5 

170 

.0 

42.84 

4.0 

80 

.0 

49, 

.11 

9, 

.0 

180 

.0 

42.45 

4.5 

90 

.0 

46 

.55 

9 

.5 

190 

.0 

42.00 

5.0 

100 

.0 

47 

.09 

10 

.0 

200 

.0 

41.62 


Case 2 

Deformation tensor = 
2.00 .00 .00 

.00 2.00 .00 

.00 .00 .25 


HRMNC 

FREQ 

PSD 

HRMMC 

FREQ 

PSO 

.5 

10 

.0 

12. 

.34 

5. 

.5 

110. 

.0 

39.34 

1.0 

20 

.0 

62. 

.08 

6 

.0 

120 

.0 

43.53 

1.5 

30 

.0 

28, 

. 14 

6 

.5 

130 

.0 

40.56 

2.0 

40 

.0 

55, 

.32 

7 

.0 

140 

.0 

41.88 

2.5 

50 

.0 

33, 

. 18 

7 

.5 

150 

.0 

41.56 

3.0 

60 

.0 

51, 

.06 

8 

.0 

160 

.0 

LO 

<s> 

<3- 

3.5 

70 

.0 

35, 

.89 

8 

.5 

170 

.0 

42.34 

4.0 

80 

.0 

47 

.95 

9 

.0 

180 

.0 

39.26 

4.5 

90 

.0 

37 

.82 

9 

.5 

190 

.0 

42.91 

5.0 

100 

.0 

45 

.51 

10 

.0 

200 

.0 

38.29 
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